knitr::opts_chunk$set(echo = TRUE)

cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")


library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(beepr)
library(DHARMa)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() beepr::beep(sound = 5)
)


####################### Model parameters #######################

iterations = 5000
warmup = 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = as.character(iterations)
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'ma[1]', 
  'cosy',
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'Daily received persuasion target -> target', 
    'Daily received persuasion target -> agent', 
    'Daily received pressure target -> target', 
    'Daily received pressure target -> agent', 
    'Daily received pushing target -> target', 
    'Daily received pushing target -> agent', 
    'Day', 
    'Daily weartime',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'Mean received persuasion target -> target',
    'Mean received persuasion target -> agent',
    'Mean received pressure target -> target',
    'Mean received pressure target -> agent',
    'Mean received pushing target -> target',
    'Mean received pushing target -> agent',
    'Mean weartime'
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Daily received persuasion target -> target)', 
  'sd(Daily received persuasion target -> agent)', 
  'sd(Daily received pressure target -> target)', 
  'sd(Daily received pressure target -> agent)', 
  'sd(Daily received pushing target -> target)', 
  'sd(Daily received pushing target -> agent)', 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]',
  'ma[1]', 
  'cosy',
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'ma[1]', 
  'cosy',
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,9),
  "Between-Person Effects" = c(10,16),
  "Random Effects" = c(17, 23), 
  "Additional Parameters" = c(24,30)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,9+5),
  "Between-Person Effects" = c(10+5,16+5),
  "Random Effects" = c(17+5, 23+5), 
  "Additional Parameters" = c(24+5,30+6)
  )

HURDLE MODELS

# For indistinguishable Dyads
model_rows_fixed_hu <- c(
    'Intercept', 
    'hu_Intercept',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    
    # HURDLE MODEL
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'hu_persuasion_self_cw', 
    'hu_persuasion_partner_cw', 
    'hu_pressure_self_cw', 
    'hu_pressure_partner_cw', 
    'hu_pushing_self_cw', 
    'hu_pushing_partner_cw', 
    'hu_day', 
    'hu_weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'hu_persuasion_self_cb',
    'hu_persuasion_partner_cb',
    'hu_pressure_self_cb',
    'hu_pressure_partner_cb',
    'hu_pushing_self_cb',
    'hu_pushing_partner_cb',
    'hu_weartime_self_cb'
  )

model_rows_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(hu_Intercept)',
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # HURDLE
  'sd(hu_persuasion_self_cw)',
  'sd(hu_persuasion_partner_cw)',
  'sd(hu_pressure_self_cw)',
  'sd(hu_pressure_partner_cw)',
  'sd(hu_pushing_self_cw)',
  'sd(hu_pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'ma[1]', 
  'cosy',
  'nu',
  'shape',
  'sderr',
  'sigma'
)
# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
    "Intercept", 
    "Hurdle Intercept",
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime",
    
    # HURDLE
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Hu Daily persuasion experienced", 
    "Hu Daily persuasion utilized (partner's view)", # OR partner received
    "Hu Daily pressure experienced", 
    "Hu Daily pressure utilized (partner's view)", 
    "Hu Daily pushing experienced", 
    "Hu Daily pushing utilized (partner's view)", 
    "Hu Day", 
    "Hu Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Hu Mean persuasion experienced", 
    "Hu Mean persuasion utilized (partner's view)", 
    "Hu Mean pressure experienced", 
    "Hu Mean pressure utilized (partner's view)", 
    "Hu Mean pushing experienced", 
    "Hu Mean pushing utilized (partner's view)", 
    "Hu Mean weartime"
  )


model_rownames_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Hurdle Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  
  # Hurdle
  "sd(Hu Daily persuasion experienced)", 
  "sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Hu Daily pressure experienced)", 
  "sd(Hu Daily pressure utilized (partner's view))", 
  "sd(Hu Daily pushing experienced)", 
  "sd(Hu Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'ma[1]', 
  'cosy',
  'nu',
  'shape',
  'sderr',
  'sigma'
)
rows_to_pack_hu <- list(
  "Conditional Within-Person Effects" = c(3,10),
  "Conditional Between-Person Effects" = c(11,17),
  
  "Hurdle Within-Person Effects" = c(18,25),
  "Hurdle Between-Person Effects" = c(26,32),
  
  "Random Effects" = c(33, 46), 
  "Additional Parameters" = c(47,53)
  )

Subjective MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Hurdle Poisson Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 1.5)", class = "b")
  , brms::set_prior("normal(0, 2.5", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(6, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  #, brms::set_prior("normal(10, 10", class = "shape")
  
  #, brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  #, brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  #, autocor_prior
)

brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = hurdle_poisson()
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  #family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  family = brms::hurdle_poisson(),
  control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_poisson_NOAR", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(pa_sub, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 12000 of 12000 iterations saturated the maximum tree depth of 10 (100%).
## Try increasing 'max_treedepth' to avoid saturation.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 322 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 12000 by 3736 log-likelihood matrix.
## 
##          Estimate     SE
## elpd_loo -32913.6 1735.6
## p_loo      4426.8  306.9
## looic     65827.1 3471.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3414  91.4%   142     
##    (0.7, 1]   (bad)       156   4.2%   <NA>    
##    (1, Inf)   (very bad)  166   4.4%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 303, observations = 3736, p-value <
## 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000267666 0.002408994
## sample estimates:
## outlier frequency (expected: 0.00133029978586724 ) 
##                                         0.08110278
#shinystan::launch_shinystan(pa_sub)

In this instance, the warning about max treedepth is a false positive. We have set treedepth to 15, and when we check with shinystan, we see that treedepth is continuously between 10 and 14.

summarize_brms(
  pa_sub, 
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_hu
)
## Warning in summarize_brms(pa_sub, model_rows_fixed = model_rows_fixed_hu, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 62.22* 55.52 69.76 1.003 1171.54 2180.69 1
Hurdle Intercept 1.20 0.87 1.66 1.005 2251.72 4451.44 0.863
Conditional Within-Person Effects
Daily persuasion experienced 0.98 0.90 1.06 1.003 820.00 1778.57 0.713
Daily persuasion utilized (partner’s view) 1.03 0.97 1.10 1.003 1094.97 2445.86 0.837
Daily pressure experienced 4.46 0.63 28.97 1.001 2038.12 3516.07 0.936
Daily pressure utilized (partner’s view) 0.88 0.72 1.07 1.001 2327.95 4209.00 0.903
Daily pushing experienced 1.05 0.89 1.27 1.002 1914.58 3765.31 0.726
Daily pushing utilized (partner’s view) 0.94 0.85 1.05 1.001 2511.46 4270.50 0.849
Day 0.92* 0.90 0.95 1.000 18563.45 9338.10 1
Daily weartime NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.25 0.94 1.69 1.001 1050.07 2075.26 0.936
Mean persuasion utilized (partner’s view) 1.01 0.75 1.37 1.001 1044.55 2129.30 0.518
Mean pressure experienced 1.04 0.78 1.39 1.001 1229.87 2676.07 0.609
Mean pressure utilized (partner’s view) 0.75* 0.56 1.00 1.001 1246.71 2682.99 0.976
Mean pushing experienced 1.14 0.76 1.71 1.003 1234.69 2416.32 0.744
Mean pushing utilized (partner’s view) 1.21 0.80 1.82 1.003 1213.71 2345.74 0.821
Mean weartime NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 0.65* 0.57 0.73 1.001 8995.78 7874.13 1
Hu Daily persuasion utilized (partner’s view) 0.75* 0.67 0.84 1.001 8190.82 7930.91 1
Hu Daily pressure experienced 1.23 0.88 1.72 1.001 9605.40 8374.53 0.897
Hu Daily pressure utilized (partner’s view) 0.66* 0.42 0.95 1.000 8978.67 6552.79 0.989
Hu Daily pushing experienced 0.58* 0.40 0.78 1.000 7451.81 7612.31 1
Hu Daily pushing utilized (partner’s view) 0.54* 0.41 0.69 1.000 7785.52 7787.54 1
Hu Day 1.09 0.84 1.42 1.000 18550.37 9136.53 0.742
Hu Daily weartime NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 0.83 0.37 1.82 1.003 2244.54 3914.95 0.682
Hu Mean persuasion utilized (partner’s view) 0.84 0.39 1.86 1.003 2220.38 3867.94 0.674
Hu Mean pressure experienced 3.32* 1.40 8.09 1.002 3571.38 5985.08 0.996
Hu Mean pressure utilized (partner’s view) 1.79 0.75 4.23 1.003 3589.16 6384.14 0.909
Hu Mean pushing experienced 0.36 0.11 1.12 1.000 3642.10 5636.41 0.964
Hu Mean pushing utilized (partner’s view) 0.34 0.11 1.06 1.000 3754.49 5708.99 0.967
Hu Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.35 0.27 0.44 1.00 2343.80 3699.48 NA
sd(Hurdle Intercept) 0.90 0.69 1.17 1.00 3156.04 5527.48 NA
sd(Daily persuasion experienced) 0.25 0.20 0.32 1.00 2676.95 4342.53 NA
sd(Daily persuasion utilized (partner’s view)) 0.19 0.15 0.25 1.00 2810.89 5109.44 NA
sd(Daily pressure experienced) 6.41 5.07 8.07 1.00 3848.63 6466.47 NA
sd(Daily pressure utilized (partner’s view)) 0.49 0.34 0.71 1.00 3944.61 6425.50 NA
sd(Daily pushing experienced) 0.51 0.33 0.76 1.00 2362.12 5293.27 NA
sd(Daily pushing utilized (partner’s view)) 0.31 0.23 0.40 1.00 4245.87 6674.84 NA
sd(Hu Daily persuasion experienced) 0.18 0.02 0.34 1.00 3802.82 2955.97 NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.17 0.02 0.33 1.00 3017.74 2932.29 NA
sd(Hu Daily pressure experienced) 0.30 0.01 0.86 1.00 4765.88 5838.00 NA
sd(Hu Daily pressure utilized (partner’s view)) 0.35 0.01 1.02 1.00 4222.93 5749.34 NA
sd(Hu Daily pushing experienced) 0.65 0.33 1.09 1.00 5623.85 7121.92 NA
sd(Hu Daily pushing utilized (partner’s view)) 0.32 0.04 0.65 1.00 4315.72 2830.48 NA
Additional Parameters
ar[1] NA NA NA NA NA NA NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr NA NA NA NA NA NA NA
sigma NA NA NA NA NA NA NA 
mcmc_plot(pa_sub,
          variable = c(
            'b_persuasion_self_cw',
            'b_persuasion_partner_cw',
            'b_pressure_self_cw',
            'b_pressure_partner_cw',
            'b_pushing_self_cw',
            'b_pushing_partner_cw'
          ),
          #regex = TRUE,
          type = 'areas',
          prob = 0.95)

mcmc_plot(pa_sub,
          variable = c(
            'b_hu_persuasion_self_cw',
            'b_hu_persuasion_partner_cw',
            'b_hu_pressure_self_cw',
            'b_hu_pressure_partner_cw',
            'b_hu_pushing_self_cw',
            'b_hu_pushing_partner_cw'
          ),
          #regex = TRUE,
          type = 'areas',
          prob = 0.95)

bayesfac <- bayestestR::bayesfactor(
  pa_sub,
  effects = "fixed"
)
## Sampling priors, please wait...
## Warning: Bayes factors might not be precise.
##   For precise Bayes factors, sampling at least 40,000 posterior samples is
##   recommended.
bayesfac
Parameter log_BF Effects Component
16 b_Intercept 189.7873193 fixed conditional
3 b_hu_Intercept -4.4904516 fixed conditional
20 b_persuasion_self_cw -3.3856527 fixed conditional
18 b_persuasion_partner_cw -3.3383318 fixed conditional
24 b_pressure_self_cw 0.6984805 fixed conditional
22 b_pressure_partner_cw -1.7927693 fixed conditional
28 b_pushing_self_cw -2.6773216 fixed conditional
26 b_pushing_partner_cw -2.7365966 fixed conditional
19 b_persuasion_self_cb -1.1528969 fixed conditional
17 b_persuasion_partner_cb -2.3378604 fixed conditional
23 b_pressure_self_cb -2.3059791 fixed conditional
21 b_pressure_partner_cb -0.2889599 fixed conditional
27 b_pushing_self_cb -1.7806387 fixed conditional
25 b_pushing_partner_cb -1.5862782 fixed conditional
1 b_day 9.8115848 fixed conditional
7 b_hu_persuasion_self_cw 11.0987333 fixed conditional
5 b_hu_persuasion_partner_cw 6.6446016 fixed conditional
11 b_hu_pressure_self_cw -1.8440165 fixed conditional
9 b_hu_pressure_partner_cw 0.0332200 fixed conditional
15 b_hu_pushing_self_cw 2.8248046 fixed conditional
13 b_hu_pushing_partner_cw 7.5653610 fixed conditional
6 b_hu_persuasion_self_cb -1.6791465 fixed conditional
4 b_hu_persuasion_partner_cb -1.6648689 fixed conditional
10 b_hu_pressure_self_cb 1.9650910 fixed conditional
8 b_hu_pressure_partner_cb -0.8626538 fixed conditional
14 b_hu_pushing_self_cb 0.1645216 fixed conditional
12 b_hu_pushing_partner_cb 0.3081348 fixed conditional
2 b_hu_day -2.7350012 fixed conditional

Additionally, as a sensitivity analysis, we estimate the two part of the models separately.

Two separate models

Modelling 0/1

The zero vs. 1 modelling part also has high pareto-k values, but reaches the same conslucsions as the hurdle model. We tried further simplifying by removing the residual AR1 correlation structure, which led to a model with good pareto-k values, still arriving at the same conslusion as the original hurdle model:

df_double$pa_sub_zero <- as.factor(ifelse(df_double$pa_sub > 0, 1, 0))

formula <- bf(
  pa_sub_zero ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
    set_prior("normal(0, 1.5)", class = "b")
  , brms::set_prior("normal(0, 5)", class = "Intercept") # for non-zero PA

  , brms::set_prior("normal(0, 1)", class = "sd", group = "coupleID")

  #, brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  #, brms::set_prior("cauchy(0, 1.5)", class = "sderr", lb = 0)
  
  #, autocor_prior
)

brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = bernoulli()
  )
## Warning: Rows containing NAs were excluded from the model.
pa_sub_zero_model <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  #family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_zero_part_NOAR", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(pa_sub_zero_model, log_pp_check = FALSE)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 2 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 12000 by 3736 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -2149.7 27.8
## p_loo        93.4  4.1
## looic      4299.5 55.6
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.3]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3734  99.9%   355     
##    (0.7, 1]   (bad)         2   0.1%   <NA>    
##    (1, Inf)   (very bad)    0   0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_sub_zero_model, integer = TRUE, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 0, observations = 3736, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0 0
## sample estimates:
## outlier frequency (expected: 0 ) 
##                                0
summarize_brms(
  pa_sub_zero_model, 
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_hu
)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 0.84 0.61 1.15 1.001 1841.97 3973.90 0.863
Hurdle Intercept NA NA NA NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced 1.53* 1.36 1.75 1.001 8062.00 8582.55 1
Daily persuasion utilized (partner’s view) 1.33* 1.19 1.50 1.000 9007.75 8833.49 1
Daily pressure experienced 0.82 0.59 1.11 1.000 10666.20 8323.82 0.906
Daily pressure utilized (partner’s view) 1.48* 1.05 2.25 1.000 10016.51 6740.30 0.987
Daily pushing experienced 1.72* 1.28 2.43 1.000 7256.09 7373.66 1
Daily pushing utilized (partner’s view) 1.83* 1.46 2.40 1.001 8929.39 8028.12 1
Day 0.92 0.71 1.19 1.000 18939.10 8954.86 0.739
Daily weartime NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.17 0.56 2.48 1.002 1891.71 4153.21 0.667
Mean persuasion utilized (partner’s view) 1.17 0.57 2.44 1.002 1970.51 4133.98 0.666
Mean pressure experienced 0.33* 0.15 0.75 1.001 2989.43 5922.42 0.996
Mean pressure utilized (partner’s view) 0.59 0.26 1.34 1.002 3409.72 5951.96 0.9
Mean pushing experienced 2.45 0.84 7.10 1.001 3185.19 5355.11 0.953
Mean pushing utilized (partner’s view) 2.56 0.88 7.30 1.001 3160.95 5396.08 0.958
Mean weartime NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA
Hu Daily pressure experienced NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA
Hu Daily pushing experienced NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA
Hu Day NA NA NA NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA
Hu Mean pressure experienced NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA
Hu Mean pushing experienced NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.89 0.68 1.14 1.00 3199.26 5243.78 NA
sd(Hurdle Intercept) NA NA NA NA NA NA NA
sd(Daily persuasion experienced) 0.18 0.02 0.34 1.00 3164.76 2298.66 NA
sd(Daily persuasion utilized (partner’s view)) 0.17 0.02 0.33 1.00 3506.22 3289.23 NA
sd(Daily pressure experienced) 0.28 0.01 0.78 1.00 4188.59 5081.55 NA
sd(Daily pressure utilized (partner’s view)) 0.31 0.01 0.90 1.00 3962.03 5276.36 NA
sd(Daily pushing experienced) 0.62 0.32 1.01 1.00 5746.41 6689.15 NA
sd(Daily pushing utilized (partner’s view)) 0.31 0.05 0.63 1.00 4373.04 3371.61 NA
sd(Hu Daily persuasion experienced) NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA
Additional Parameters
ar[1] NA NA NA NA NA NA NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr NA NA NA NA NA NA NA
sigma NA NA NA NA NA NA NA 
mcmc_plot(pa_sub_zero_model,
          variable = c(
            'b_persuasion_self_cw',
            'b_persuasion_partner_cw',
            'b_pressure_self_cw',
            'b_pressure_partner_cw',
            'b_pushing_self_cw',
            'b_pushing_partner_cw'
          ),
          #regex = TRUE,
          type = 'areas',
          prob = 0.95)

Modelling active Days

# Only include active days
df_double$pa_sub_non_zero <- ifelse(df_double$pa_sub > 0, df_double$pa_sub, NA)

formula <- bf(
  log(pa_sub_non_zero) ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  ##, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA

  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  , brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  #, brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  #, autocor_prior
)

brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = gaussian()
  )
## Warning: Rows containing NAs were excluded from the model.
pa_sub_active_model <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_active_part_NOAR", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(pa_sub_active_model, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 12000 by 1672 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo  -1785.5 35.3
## p_loo        85.9  4.1
## looic      3571.1 70.6
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.0]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_sub_active_model, integer = TRUE, outliers_type = 'bootstrap')
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 8, observations = 1672, p-value = 0.26
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0005980861 0.0068929426
## sample estimates:
## outlier frequency (expected: 0.00320574162679426 ) 
##                                        0.004784689
summarize_brms(
  pa_sub_active_model, 
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_hu
)
## Warning in summarize_brms(pa_sub_active_model, model_rows_fixed =
## model_rows_fixed_hu, : Coefficients were exponentiated. Double check if this
## was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 47.85* 42.22 54.17 1.001 4041.60 6381.47 1
Hurdle Intercept NA NA NA NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced 1.03 0.97 1.08 1.000 6597.52 8346.37 0.833
Daily persuasion utilized (partner’s view) 1.03 0.99 1.08 1.000 8955.61 9602.22 0.904
Daily pressure experienced 0.89* 0.80 0.99 1.000 14496.32 8806.27 0.987
Daily pressure utilized (partner’s view) 0.94 0.85 1.03 1.000 13817.24 8187.83 0.908
Daily pushing experienced 1.03 0.96 1.10 1.000 11709.14 9409.07 0.775
Daily pushing utilized (partner’s view) 0.99 0.93 1.05 1.000 11995.74 10105.19 0.633
Day 1.01 0.89 1.14 1.000 20035.93 9281.40 0.545
Daily weartime NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.00 0.73 1.37 1.001 3231.81 4989.35 0.51
Mean persuasion utilized (partner’s view) 0.98 0.72 1.33 1.001 3286.72 5489.83 0.56
Mean pressure experienced 1.15 0.80 1.64 1.000 5002.27 7220.72 0.778
Mean pressure utilized (partner’s view) 0.89 0.61 1.28 1.000 5018.40 7058.04 0.732
Mean pushing experienced 1.34 0.84 2.11 1.000 4624.61 7355.15 0.896
Mean pushing utilized (partner’s view) 1.42 0.89 2.26 1.000 4706.03 7431.81 0.929
Mean weartime NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA
Hu Daily pressure experienced NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA
Hu Daily pushing experienced NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA
Hu Day NA NA NA NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA
Hu Mean pressure experienced NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA
Hu Mean pushing experienced NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.32 0.25 0.42 1.00 3621.92 6455.57 NA
sd(Hurdle Intercept) NA NA NA NA NA NA NA
sd(Daily persuasion experienced) 0.12 0.08 0.17 1.00 7194.08 9303.03 NA
sd(Daily persuasion utilized (partner’s view)) 0.09 0.05 0.13 1.00 8556.16 8646.99 NA
sd(Daily pressure experienced) 0.08 0.00 0.23 1.00 5748.83 6157.17 NA
sd(Daily pressure utilized (partner’s view)) 0.07 0.00 0.19 1.00 6534.51 6674.27 NA
sd(Daily pushing experienced) 0.11 0.04 0.19 1.00 5752.73 4564.11 NA
sd(Daily pushing utilized (partner’s view)) 0.09 0.02 0.17 1.00 5084.36 3390.77 NA
sd(Hu Daily persuasion experienced) NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA
Additional Parameters
ar[1] NA NA NA NA NA NA NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr NA NA NA NA NA NA NA
sigma 0.68 0.66 0.71 1.00 19220.03 8549.91 NA 
mcmc_plot(pa_sub_active_model,
          variable = c(
            'b_persuasion_self_cw',
            'b_persuasion_partner_cw',
            'b_pressure_self_cw',
            'b_pressure_partner_cw',
            'b_pushing_self_cw',
            'b_pushing_partner_cw'
          ),
          #regex = TRUE,
          type = 'areas',
          prob = 0.95)

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

log gaussian

formula <- bf(
  log(pa_obj) ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA

  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  , brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  #, brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  #, autocor_prior
)

brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = gaussian()
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_NOAR", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(pa_obj_log, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 12000 by 3337 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2936.5  56.3
## p_loo        92.1   4.5
## looic      5872.9 112.6
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 26, observations = 3337, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.001341025 0.004952053
## sample estimates:
## outlier frequency (expected: 0.00295474977524723 ) 
##                                        0.007791429
summarize_brms(
  pa_obj_log, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Warning in summarize_brms(pa_obj_log, model_rows_fixed = model_rows_fixed, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 117.41* 105.48 130.41 1.002 1936.20 4337.49 1
Within-Person Effects
Daily persuasion experienced 1.03 1.00 1.06 1.000 8443.32 9440.38 0.966
Daily persuasion utilized (partner’s view) 1.02 0.99 1.05 1.001 10110.26 9523.00 0.888
Daily pressure experienced 0.94 0.88 1.01 1.000 13123.29 8882.26 0.96
Daily pressure utilized (partner’s view) 0.98 0.92 1.05 1.000 14150.05 9457.65 0.714
Daily pushing experienced 1.03 0.98 1.08 1.000 11136.01 8673.71 0.9
Daily pushing utilized (partner’s view) 1.02 0.97 1.06 1.000 15288.78 9998.46 0.771
Day 0.97 0.91 1.04 1.000 20771.69 8890.89 0.785
Daily weartime 1.00* 1.00 1.00 1.000 12258.92 8147.17 1
Between-Person Effects
Mean persuasion experienced 1.10 0.82 1.46 1.003 1635.21 3066.39 0.748
Mean persuasion utilized (partner’s view) 0.98 0.73 1.30 1.003 1664.34 3415.83 0.559
Mean pressure experienced 0.98 0.73 1.31 1.003 2233.86 4590.15 0.56
Mean pressure utilized (partner’s view) 0.97 0.73 1.28 1.002 2141.33 4170.53 0.594
Mean pushing experienced 0.97 0.65 1.45 1.002 2615.53 4583.92 0.553
Mean pushing utilized (partner’s view) 1.25 0.83 1.85 1.002 2546.44 4664.09 0.859
Mean weartime 1.00 1.00 1.00 1.000 16930.39 10509.93 0.909
Random Effects
sd(Intercept) 0.31 0.24 0.40 1.00 2373.00 4540.45 NA
sd(Daily persuasion experienced) 0.05 0.02 0.08 1.00 6108.14 4438.88 NA
sd(Daily persuasion utilized (partner’s view)) 0.06 0.03 0.09 1.00 6116.29 5901.98 NA
sd(Daily pressure experienced) 0.05 0.00 0.14 1.00 5104.76 5590.89 NA
sd(Daily pressure utilized (partner’s view)) 0.04 0.00 0.12 1.00 6898.74 5731.73 NA
sd(Daily pushing experienced) 0.07 0.01 0.15 1.00 2745.29 3354.31 NA
sd(Daily pushing utilized (partner’s view)) 0.04 0.00 0.10 1.00 4176.37 5266.68 NA
Additional Parameters
ar[1] NA NA NA NA NA NA NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr NA NA NA NA NA NA NA
sigma 0.57 0.56 0.59 1.00 19512.21 8255.47 NA 
mcmc_plot(pa_obj_log,
          variable = c(
            'b_persuasion_self_cw',
            'b_persuasion_partner_cw',
            'b_pressure_self_cw',
            'b_pressure_partner_cw',
            'b_pushing_self_cw',
            'b_pushing_partner_cw'
          ),
          #regex = TRUE,
          type = 'areas',
          prob = 0.95)

Affect

range(df_double$aff, na.rm = T) 
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6) # range of the outcome scale
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  , brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  #, brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  #, autocor_prior
)

brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = gaussian()
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_NOAR", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(mood_gauss, log_pp_check = FALSE)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 12000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -5185.9  59.2
## p_loo        75.2   3.3
## looic     10371.7 118.5
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.1]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(mood_gauss, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 31, observations = 3736, p-value =
## 9.752e-11
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.00564461 0.01175733
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.008297645
summarize_brms(
  mood_gauss, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 3.70* 3.48 3.91 1.007 1299.69 2626.71 1
Within-Person Effects
Daily persuasion experienced 0.00 -0.04 0.05 1.000 10096.58 8360.31 0.553
Daily persuasion utilized (partner’s view) 0.02 -0.02 0.07 1.001 8776.31 8901.48 0.829
Daily pressure experienced -0.04 -0.14 0.07 1.000 11418.22 8197.57 0.767
Daily pressure utilized (partner’s view) -0.02 -0.14 0.08 1.000 10703.53 8437.61 0.661
Daily pushing experienced 0.02 -0.04 0.09 1.000 11578.47 8985.61 0.768
Daily pushing utilized (partner’s view) 0.08* 0.01 0.14 1.000 10076.04 8287.53 0.984
Day 0.26* 0.15 0.37 1.001 17305.49 8858.10 1
Daily weartime NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.33 -0.21 0.90 1.003 1324.40 2167.41 0.889
Mean persuasion utilized (partner’s view) 0.22 -0.32 0.79 1.003 1322.88 2252.89 0.79
Mean pressure experienced -0.31 -0.86 0.23 1.003 1481.77 2769.72 0.874
Mean pressure utilized (partner’s view) -0.30 -0.84 0.24 1.003 1476.58 2758.08 0.871
Mean pushing experienced 0.22 -0.54 0.99 1.002 1997.17 3840.31 0.714
Mean pushing utilized (partner’s view) 0.35 -0.39 1.12 1.002 2005.43 3822.38 0.825
Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.60 0.47 0.78 1.00 2148.52 4074.04 NA
sd(Daily persuasion experienced) 0.04 0.00 0.10 1.00 3348.96 5109.75 NA
sd(Daily persuasion utilized (partner’s view)) 0.08 0.01 0.13 1.00 3084.49 2127.56 NA
sd(Daily pressure experienced) 0.08 0.00 0.24 1.00 4964.63 5526.49 NA
sd(Daily pressure utilized (partner’s view)) 0.09 0.00 0.26 1.00 4911.35 5368.11 NA
sd(Daily pushing experienced) 0.06 0.00 0.14 1.00 3656.74 3408.48 NA
sd(Daily pushing utilized (partner’s view)) 0.07 0.00 0.17 1.00 4184.55 4073.31 NA
Additional Parameters
ar[1] NA NA NA NA NA NA NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr NA NA NA NA NA NA NA
sigma 0.96 0.94 0.98 1.00 18080.51 9327.96 NA 
mcmc_plot(mood_gauss,
          variable = c(
            'b_persuasion_self_cw',
            'b_persuasion_partner_cw',
            'b_pressure_self_cw',
            'b_pressure_partner_cw',
            'b_pushing_self_cw',
            'b_pushing_partner_cw'
          ),
          #regex = TRUE,
          type = 'areas',
          prob = 0.95)

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  #, brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) # range of the outcome scale
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  #, brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  #, brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  #, autocor_prior
)


brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = cumulative() # HURDLE_CUMULATIVE
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal_NOARNOAR_", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(reactance_ordinal)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 6 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -681.9 31.9
## p_loo        73.4  5.4
## looic      1363.8 63.8
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.7]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     750   99.2%   445     
##    (0.7, 1]   (bad)        6    0.8%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 2, observations = 756, p-value = 0.08
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.002645503
## sample estimates:
## outlier frequency (expected: 0.00041005291005291 ) 
##                                        0.002645503
summarize_brms(
  reactance_ordinal, 
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercepts
Intercept NA NA NA NA NA NA NA
Intercept[1] 3.85* 2.33 6.45 1.000 8597.51 9211.52 1
Intercept[2] 8.35* 4.95 14.45 1.000 8834.13 8761.87 1
Intercept[3] 23.24* 13.13 42.31 1.000 9373.11 9135.25 1
Intercept[4] 101.58* 52.10 209.10 1.000 10595.59 9691.96 1
Intercept[5] 3488.39* 1077.21 13336.21 1.001 13289.90 8614.35 1
Within-Person Effects
Daily persuasion experienced 0.85* 0.71 0.99 1.000 9896.79 7775.98 0.98
Daily persuasion utilized (partner’s view) 1.03 0.84 1.24 1.000 9229.64 7954.31 0.607
Daily pressure experienced 1.84* 1.18 2.68 1.001 5608.70 6574.26 0.994
Daily pressure utilized (partner’s view) 1.22 0.71 2.06 1.001 7420.47 6712.38 0.808
Daily pushing experienced 1.17 0.97 1.43 1.001 7757.54 7485.59 0.946
Daily pushing utilized (partner’s view) 0.91 0.71 1.17 1.000 9792.98 8248.45 0.77
Day 1.47 0.75 2.89 1.000 13856.59 9556.08 0.871
Daily weartime NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.12 0.40 3.11 1.000 4004.78 6205.61 0.586
Mean persuasion utilized (partner’s view) 1.38 0.45 4.31 1.000 4278.63 6604.64 0.711
Mean pressure experienced 3.51* 1.18 10.71 1.000 4792.44 6629.78 0.99
Mean pressure utilized (partner’s view) 1.17 0.37 3.66 1.000 4762.77 6973.89 0.612
Mean pushing experienced 1.23 0.28 5.62 1.000 5266.91 7722.24 0.605
Mean pushing utilized (partner’s view) 0.11* 0.02 0.64 1.000 7169.66 7894.75 0.992
Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.82 0.47 1.25 1.00 4256.66 7026.55 NA
sd(Daily persuasion experienced) 0.18 0.01 0.43 1.00 1818.16 4469.30 NA
sd(Daily persuasion utilized (partner’s view)) 0.22 0.01 0.51 1.00 3273.96 4646.21 NA
sd(Daily pressure experienced) 0.56 0.09 1.14 1.00 2689.26 2386.75 NA
sd(Daily pressure utilized (partner’s view)) 0.50 0.02 1.55 1.00 2929.20 4951.16 NA
sd(Daily pushing experienced) 0.22 0.01 0.50 1.00 3370.38 4130.06 NA
sd(Daily pushing utilized (partner’s view)) 0.18 0.01 0.57 1.00 4830.97 5284.06 NA
Additional Parameters
ar[1] NA NA NA NA NA NA NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr NA NA NA NA NA NA NA
sigma NA NA NA NA NA NA NA
disc 1.00 1.00 1.00 NA NA NA NA 
mcmc_plot(reactance_ordinal,
          variable = c(
            'b_persuasion_self_cw',
            'b_persuasion_partner_cw',
            'b_pressure_self_cw',
            'b_pressure_partner_cw',
            'b_pushing_self_cw',
            'b_pushing_partner_cw'
          ),
          #regex = TRUE,
          type = 'areas',
          prob = 0.95)

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) # range of the outcome scale
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  #, brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  #, brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  #, autocor_prior
)


brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = bernoulli()
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_NOARNOAR_", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(is_reactance)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 55 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -363.1 16.0
## p_loo        80.1  6.0
## looic       726.3 32.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.5]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     701   92.7%   429     
##    (0.7, 1]   (bad)       49    6.5%   <NA>    
##    (1, Inf)   (very bad)   6    0.8%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(is_reactance, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0000334886 0.0073476538
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.001322751
summarize_brms(
  is_reactance, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 0.29* 0.16 0.50 1.000 13412.80 10197.69 1
Within-Person Effects
Daily persuasion experienced 0.84 0.69 1.01 1.001 12653.33 8961.78 0.966
Daily persuasion utilized (partner’s view) 1.13 0.85 1.54 1.000 9878.26 8807.43 0.789
Daily pressure experienced 2.04* 1.03 4.57 1.000 8289.60 6530.31 0.979
Daily pressure utilized (partner’s view) 1.44 0.58 4.04 1.000 8810.48 6987.21 0.799
Daily pushing experienced 1.28* 1.01 1.64 1.001 12865.46 8547.86 0.979
Daily pushing utilized (partner’s view) 0.89 0.60 1.31 1.000 14162.09 8870.67 0.735
Day 1.64 0.77 3.51 1.000 19648.18 9310.83 0.9
Daily weartime NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.99 0.60 6.80 1.000 6422.05 7842.49 0.87
Mean persuasion utilized (partner’s view) 1.90 0.52 7.15 1.000 6799.76 8735.29 0.829
Mean pressure experienced 17.91* 2.33 159.59 1.000 8975.03 8541.37 0.997
Mean pressure utilized (partner’s view) 2.27 0.24 19.14 1.000 8119.78 8948.63 0.769
Mean pushing experienced 0.83 0.12 6.30 1.000 8730.85 8605.26 0.58
Mean pushing utilized (partner’s view) 0.08* 0.01 0.66 1.000 10022.27 10004.57 0.99
Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.18 0.75 1.74 1.00 5382.54 7890.40 NA
sd(Daily persuasion experienced) 0.21 0.01 0.51 1.00 2939.60 5289.74 NA
sd(Daily persuasion utilized (partner’s view)) 0.50 0.12 0.97 1.00 4502.23 4854.72 NA
sd(Daily pressure experienced) 1.12 0.16 2.43 1.00 2982.38 3385.24 NA
sd(Daily pressure utilized (partner’s view)) 0.97 0.04 2.73 1.00 4249.53 5650.95 NA
sd(Daily pushing experienced) 0.25 0.01 0.60 1.00 4754.94 5282.30 NA
sd(Daily pushing utilized (partner’s view)) 0.31 0.01 0.95 1.00 4810.38 6592.98 NA
Additional Parameters
ar[1] NA NA NA NA NA NA NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr NA NA NA NA NA NA NA
sigma NA NA NA NA NA NA NA 
mcmc_plot(is_reactance,
          variable = c(
            'b_persuasion_self_cw',
            'b_persuasion_partner_cw',
            'b_pressure_self_cw',
            'b_pressure_partner_cw',
            'b_pushing_self_cw',
            'b_pushing_partner_cw'
          ),
          #regex = TRUE,
          type = 'areas',
          prob = 0.95)

hypothesis(is_reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0     0.47      0.39    -0.15     1.13       8.86
##   Post.Prob Star
## 1       0.9     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

if (report_ordinal) {
  model_rows_random_final <- model_rows_random_ordinal
  model_rows_fixed_final <- model_rows_fixed_ordinal
  model_rownames_fixed_final <- model_rownames_fixed_ordinal
  model_rownames_random_final <- model_rownames_random_ordinal
  rows_to_pack_final <- rows_to_pack_ordinal
} else {
  model_rows_random_final <- model_rows_random_hu
  model_rows_fixed_final <- model_rows_fixed_hu
  model_rownames_fixed_final <- model_rownames_fixed_hu
  model_rownames_random_final <- model_rownames_random_hu
  rows_to_pack_final <- rows_to_pack_hu
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  model_rows_random = model_rows_random_final,
  model_rows_fixed = model_rows_fixed_final,
  model_rownames_random = model_rownames_random_final,
  model_rownames_fixed = model_rownames_fixed_final
) 
## [1] "pa_sub"
## Warning in summarize_brms(model, short_version = TRUE, exponentiate =
## exponentiate, : Coefficients were exponentiated. Double check if this was
## intended.
## [1] "pa_obj_log"
## Warning in summarize_brms(model, short_version = TRUE, exponentiate =
## exponentiate, : Coefficients were exponentiated. Double check if this was
## intended.
## [1] "mood_gauss"
## [1] "reactance_ordinal"
## [1] "is_reactance"
# pretty printing
summary_all_models <- all_models %>%
  print_df(rows_to_pack = rows_to_pack_final) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Poisson" = 2, 
      "Device-Based MVPA Log (Gaussian)" = 2, 
      "Mood Gaussian" = 2,
      "Reactance Ordinal" = 2,
      "Reactance Dichotome" = 2)
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack_final,
  file.path("Output", "AllModels_experimental_noAR.xlsx"), 
  merge_option = 'both', 
  simplify_2nd_row = TRUE,
  colwidths = c(38, 7.2, 13.3, 7.2, 13.3,7.2, 13.3,7.2, 13.3,7.2, 13.3),
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
summary_all_models
Subjective MVPA Hurdle Poisson
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Ordinal
Reactance Dichotome
exp(Est.) pa_sub p_direction pa_sub exp(Est.) pa_obj_log p_direction pa_obj_log b mood_gauss p_direction mood_gauss OR reactance_ordinal p_direction reactance_ordinal OR is_reactance p_direction is_reactance
Intercept 62.22* 1 117.41* 1 3.70* 1 NA NA 0.29* 1
Hurdle Intercept 1.20 0.863 NA NA NA NA NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced 0.98 0.713 1.03 0.966 0.00 0.553 0.85* 0.98 0.84 0.966
Daily persuasion utilized (partner’s view) 1.03 0.837 1.02 0.888 0.02 0.829 1.03 0.607 1.13 0.789
Daily pressure experienced 4.46 0.936 0.94 0.96 -0.04 0.767 1.84* 0.994 2.04* 0.979
Daily pressure utilized (partner’s view) 0.88 0.903 0.98 0.714 -0.02 0.661 1.22 0.808 1.44 0.799
Daily pushing experienced 1.05 0.726 1.03 0.9 0.02 0.768 1.17 0.946 1.28* 0.979
Daily pushing utilized (partner’s view) 0.94 0.849 1.02 0.771 0.08* 0.984 0.91 0.77 0.89 0.735
Day 0.92* 1 0.97 0.785 0.26* 1 1.47 0.871 1.64 0.9
Daily weartime NA NA 1.00* 1 NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.25 0.936 1.10 0.748 0.33 0.889 1.12 0.586 1.99 0.87
Mean persuasion utilized (partner’s view) 1.01 0.518 0.98 0.559 0.22 0.79 1.38 0.711 1.90 0.829
Mean pressure experienced 1.04 0.609 0.98 0.56 -0.31 0.874 3.51* 0.99 17.91* 0.997
Mean pressure utilized (partner’s view) 0.75* 0.976 0.97 0.594 -0.30 0.871 1.17 0.612 2.27 0.769
Mean pushing experienced 1.14 0.744 0.97 0.553 0.22 0.714 1.23 0.605 0.83 0.58
Mean pushing utilized (partner’s view) 1.21 0.821 1.25 0.859 0.35 0.825 0.11* 0.992 0.08* 0.99
Mean weartime NA NA 1.00 0.909 NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 0.65* 1 NA NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) 0.75* 1 NA NA NA NA NA NA NA NA
Hu Daily pressure experienced 1.23 0.897 NA NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) 0.66* 0.989 NA NA NA NA NA NA NA NA
Hu Daily pushing experienced 0.58* 1 NA NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) 0.54* 1 NA NA NA NA NA NA NA NA
Hu Day 1.09 0.742 NA NA NA NA NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 0.83 0.682 NA NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) 0.84 0.674 NA NA NA NA NA NA NA NA
Hu Mean pressure experienced 3.32* 0.996 NA NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) 1.79 0.909 NA NA NA NA NA NA NA NA
Hu Mean pushing experienced 0.36 0.964 NA NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) 0.34 0.967 NA NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.35 NA 0.31 NA 0.60 NA 0.82 NA 1.18 NA
sd(Hurdle Intercept) 0.90 NA NA NA NA NA NA NA NA NA
sd(Daily persuasion experienced) 0.25 NA 0.05 NA 0.04 NA 0.18 NA 0.21 NA
sd(Daily persuasion utilized (partner’s view)) 0.19 NA 0.06 NA 0.08 NA 0.22 NA 0.50 NA
sd(Daily pressure experienced) 6.41 NA 0.05 NA 0.08 NA 0.56 NA 1.12 NA
sd(Daily pressure utilized (partner’s view)) 0.49 NA 0.04 NA 0.09 NA 0.50 NA 0.97 NA
sd(Daily pushing experienced) 0.51 NA 0.07 NA 0.06 NA 0.22 NA 0.25 NA
sd(Daily pushing utilized (partner’s view)) 0.31 NA 0.04 NA 0.07 NA 0.18 NA 0.31 NA
sd(Hu Daily persuasion experienced) 0.18 NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.17 NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) 0.30 NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) 0.35 NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) 0.65 NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) 0.32 NA NA NA NA NA NA NA NA NA
Additional Parameters
ar[1] NA NA NA NA NA NA NA NA NA NA
ma[1] NA NA NA NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA NA NA NA
sderr NA NA NA NA NA NA NA NA NA NA
sigma NA NA 0.57 NA 0.96 NA NA NA NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::report_packages()
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.26.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.21.0; Bürkner P, 2017)
  • Rcpp (version 1.0.13; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.6; Hartig F, 2022)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • R (version 4.4.1; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • knitr (version 1.48; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()